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In this paper we describe the semi-spectral linear MHD (SLiM) code which we have written to follow the interaction 
of linear waves through an inhomogeneous three-dimensional solar atmosphere. The background model allows almost 
arbitrary perturbations of density, temperature, sound speed as well as magnetic and velocity fields. We give details of 
several of the tests we have used to check the code. The code will be useful in understanding the helioseismic signatures 
of various solar features, including sunspots. 



1 Introduction 



Local helioseismology uses waves propagating through the 
solar interior to probe the inhomogeneities therein. Flows, 
temperature and magnetic field inhomogeneities all imprint 
their signatures on the waves, and helioseismology attempts 
to learn about the inhomogeneities by studying the waves 
(Duvall et al. 1993). In order to better understand the various 
signatures we have developed a numerical code for forward 
calculations which numerically propagates waves through 
various types of inhomogeneities. 

Such calculations have been performed in the past, but 
mainly in a 2-D geometry (eg Cally & Bogdan 1997), or 
with other simplifying physical assumptions (eg Bogdan et 
al. 1996; Fan, Braun & Chou 1995). Analytic solutions also 
exist for some idealised problems: sound speed perturba- 
tions (eg Jensen, & Pijpers 2003), an isolated mass flow 
(Birch & Felder 2004), or a flux tube in an otherwise uni- 
form atmosphere (Gizon, Hanasoge, & Birch 2006). These 
solutions have been very useful in revealing various effects 
and in aiding our understanding. There remains however 
a need for 3 -dimensional treatments (see Werne, Birch & 
JuUen 2004). 

This paper describes the three-dimensional physical and 
mathematical problem we have chosen to tackle, as well as 
the numerical scheme we have implemented. The code has 
been subjected to a number of tests, several of which will 
be described in detail. Our intention is to make the code 
available to the helioseismology community as part of the 
HELAS project. 
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2 Statement of the problem 

2.1 The nonlinear equations 

As our starting point we consider a typical formulation of 
the MHD equations in cgs units. The relevant equations are 
the continuity equation, 

dtp + V ■{pU) = 0, (1) 

the compressible Navier-Stokes equation including Lorentz 
and buoyancy forces, 

pDtU = -VP + pgz + ^{V X B) X B 



4tt 



fv(v.c/). 



(2) 



(3) 



(4) 
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where Dt = dt + U ■ 'V, the induction equation 
dtB^V x{U xB)+ riAB , 
an energy equation, 
pCyDtT = -PV ■ U + KV^T 

+viscous heating -I- ohmic heating 

+radiative heating , 

an equation of state, P — P{p,T), and the solenoidality 
condition V ■ B = 0. Usual symbols are used to denote the 
physical quantities that appear in the above equations. The 
gravitational acceleration, g, is assumed to be uniform and 
constant. The diffusivities as written here are particularly 
simple, being both isotropic and uniform, and in principle 
more complicated forms might be needed in certain appli- 
cations, however as we intend to treat the waves as ideal 
perturbations and ignore the effects of the diffusivities, the 
particular forms of these terms is not important. 

2.2 The ideal, Unear equations 

The waves are generally small perturbations, excited by the 
turbulent motions of the granules or in rare cases as a re- 
sult of impulsive events such as flares. Our code follows 
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the evolution of a such perturbation superimposed on a pro- 
vided background. The background, assumed to be time- 
independent and in equihbrium, and is specified by giving 
the density po{r), pressure Po{r), magnetic field Bo(r), 
velocity field Uo{r) and the first adiabatic index ri(r) at 
position r = {x, y, z) everywhere inside the domain. The 
background current, Jq = V x i?o, and sound speed, 
Co ~ (ri-Po/po)^^^. can be derived from the above quanti- 
ties. 

The perturbation is assumed to be sufficiently small that 
we can consider the linearised equations. We currently as- 
sume the waves are adiabatic and the wave propagation is 
ideal. This allows us to write the magnetic, velocity, pres- 
sure and density perturbations in terms of the displacement 
vector ^. Using primed quantities as the Eulerian perturba- 
tions to the background state, the master equation for ^ is 

p^d'ii = -VP' + p'gz + ^{J' X Bo + Jo X B') 

Ait 



+ {podti - p'Uo - PoU') ■ VUo 
-poUo-V{U + U'), 



(5) 



where the perturbed quantities that appear in this equation 
are themselves functions of ^ according to 



p' - -V.(poO, 

P' - clip' + ^-Vpo)-^-VPo, 

B' = Vxi^x Bo), 

J' = V X B', 

U' = dt^ + Uo-V^-^-VUo. 



(6) 
(7) 
(8) 
(9) 
(10) 



We refer the reader to Lynden-Bell and Ostriker (1967) for 
equation ( fTOl i. to Christensen-Dalsgaard (2003, page 174) 
for a treatment of the continuity equation (eq. ||6l), and to 
Moffatt (1978) for the induction equation (eq. [S]). Note that 
we have assumed that the presence of a flow and a magnetic 
field are mutually exclusive (eq. JH]). 

2.3 The boundary conditions 

The tests described in this paper involve magnetic, density 
or velocity inhomogeneities superimposed on an isotropic, 
uniform background. For these tests we have then assumed 
periodicity in all three spatial directions. However since the 
code is also intended to be used for stratified atmospheres, 
we will here describe the boundary condition we have im- 
plemented for such cases. 

We only intend to simulate small regions of the solar 
surface and so have used a box-geometry and used a spec- 
tral treatment in the horizontal directions and finite differ- 
ence for the vertical. The boundary conditions are naturally 
important. The side boundaries are simple: the box is peri- 
odic. For the upper boundary we currently use the condition 
that the Lagrangian perturbation of the vertical component 
of the stress tensor vanishes (following Cally and Bogdan 
1997): 



where 11 is the stress tensor. Furthermore we assume that the 
background flow is restricted to regions away from the top 
boundary, and that there is no purely horizontal field there, 
so that 



(12) 



In regions where Bo 7^ the condition + ^ • ^^xz = 
imphes 



B'=- 



BoxB'^ + $ ■ W{BoxBoz 



Oz 



and the condition on Il'y^ implies 



b:. 



BpyB'^ + ^ • ^{BpyBpz) 
Boz 



(13) 



(14) 



The boundary condition 11'^^ + ^ • Vn^^ = can then be 
converted into a constraint on P' at the top boundary: 



P' = ^i^OzB', - BoxB', - BoyB'y) - ^ • VPo 
-j^i-'^iBlx+Bly-Bl,). 



(15) 



Once this is known we can also deduce the density at the 
top boundary using equation (7): 



p' - cp\p' + i-vPo)-i-vpo. 



(16) 



Lastly, our assumption that the velocity vanishes near the 
boundary implies 



u' - dtS,. 



(17) 



n^, + 1 • vn,, = 0, 



(11) 



In those regions where Bo — 0, the boundary conditions 
have B' — Q but are otherwise unchanged. 

In the case where there is a magnetic field at the bound- 
ary, it can readily be seen from equation (8) that specifying 
B'^ on the boundary is equivalent to specifying a relation- 
ship between dz£,x and dz£,z- Similarly a constraint on dz£^y 
and dz£,z is given by the condition on B'y and all three are 
related by dzS,z by the condition on p' and equation (6). We 
thus have three relationships for three unknowns and have 
a well posed problem for defining 9^^: the boundary con- 
ditions, in magnetic regions, are equivalent to specifying a 
boundary condition three components of dz^,- 

In the non-magnetic there are fewer wave modes and 
the boundary condition reduces to specifying dz£,z- In this 
special case the boundary condition corresponds to that of a 
free-surface (Cally and Bogdan 1997). 

2.4 Initial Conditions 

The above equations allow us to follow the evolution of 
arbitrary perturbations ^(r,to) and dt^{r,to)- Often these 
initial conditions will take the form of an f- or p-mode 
packet, but our code is able to deal with arbitrary initial per- 
turbations. 
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2.5 Physical instabilities 

The code evolves perturbations superimposed on the back- 
ground atmosphere assuming that the evolution is linear. 
Consequently if the atmosphere is (linearly) unstable then 
exponentially growing modes exist. Under almost all cir- 
cumstances the fastest growing of these modes will even- 
tually dominate the solution. This is particularly important 
to remember when our interest is wave propagation in the 
convection zone where the atmosphere is in places supera- 
diabatically stratified. In such places the atmosphere is con- 
vectively unstable, and exponentially growing solutions are 
to be expected: this is a real property of the background, and 
will necessarily exist in any proper numerical treatment of 
the problem. 

Our suggestion, and the approach we have adopted in 
some preliminary calculations, is to preprocess the back- 
ground so that it becomes marginally stable. This removes 
the unstable modes without introducing new unphysical 
modes, for example gravity waves in the convection zone. 
Naturally the code can also be used to study the linear 
growth of ideal instabilities without modification to either 
the code or the atmosphere. 



3 Numerical Methods 

Since we have chosen periodic boundary conditions for the 
side walls, we have implemented a pseudo-spectral scheme 
in the horizontal directions. That is all horizontal derivatives 
are evaluated in Fourier space, all products are evaluated in 
physical space: we use the FFTW package (Frigo & John- 
son 2005) to move between these two spaces. Since our aim 
is to model the Sun, stratification due to gravity imphes pe- 
riodic boundaries are not appropriate for the top and bot- 
tom boundaries. We have therefore used a simple two-step 
Lax-Wendroff scheme in the vertical to evolve the horizon- 
tal Fourier modes. Care has been taken in order to ensure 
that the various derivatives are defined on the appropriate 
grids. 

Since the Lax-Wendroff scheme is reasonably efficient 
in damping waves with wavelengths comparable to the grid- 
spacing, we have implemented a hyper-diffusivity in the 
horizontal directions which scales like fc^, where k is the 
horizontal wavenumber The short wavelength modes in the 
horizontal directions are thus damped in a similar way to 
that which occurs in the vertical without significantly af- 
fecting those modes we are interested in. 



4 Comparison with exact solutions 

Here we provide some examples of test calculations used to 
verify the code. For these tests, the background is uniform 
away from the inhomogeneity and the initial condition is a 
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Fig. 1 Scattered pressure perturbations (snapshot) from a 
vertical 1 kG flux tube with axis x = y = 0. The radius of 
the tube is 2 Mm. The bottom panel shows the exact result 
and the top panel the numerical result. 

simple acoustic wave packet propagating in the +x direc- 
tion (to the right): 

/>oo 

|(r, to) = -x A{k) sm[k{x - xq)] dk, (18) 
Jo 

and 

poo 

dt^{r,to) = X / A{k) cok cos[k{x — xo)]dk, (19) 
Jo 

where r — (x, y. z) and A{k) is a Gaussian envelope cen- 
tred on the wavenumber corresponding to angular frequency 
3 mHz and with standard deviation 1 mHz. Att = to the 
wave packet is centred at the left edge of the computational 
domain with xq = —12 Mm. 

4.1 Magnetic and pressure perturbations 

Our first test reproduces the analytic result of Wilson (1980) 
and Gizon et al. (2006) where the problem is fuUy de- 
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Fig. 2 A cut at y = through the scattered pressure per- 
turbations shown in Figure [T] (vertical 1 kG flux tube). The 
dashed line shows the exact solution and the solid line the 
numerical. The units are such that the incoming pressure 
perturbations have an amplitude of 1 . 

scribed. In brief, the background state consists of a cylin- 
drical flux tube in pressure and thermal equilibrium. The 
external plasma is to be uniform and at rest, with density 
5 X 10"^ g/cm^, sound speed 11 x 10'' cm/s, and with 
Fi = 5/3. The flux tube, which has a radius of 2 Mm in our 
tests, has uniform properties. Its magnetic field is uniform, 
and we have considered different field strengths -Btubc from 
1 kG (relatively weak) to 3 kG (where the flux tube is 99% 
evacuated). The magnetic field and the condition of pressure 
equilibrium imply that the pressure inside the tube is 

Pin = Pcxt - 5tlbc/(87r). (20) 

The assumed thermal equilibrium requires that the tempera- 
ture inside the tube equals that outside the tube and that the 
density inside the tube be reduced. 

The problem thus has a discontinuous jump across the 
flux tube. This jump is important for obtaining the analytic 
solution but is somewhat problematic for a numerical treat- 
ment. Fortunately there is no reason for expecting the dis- 
continuity to be singular, and we can therefore proceed by 
considering the discontinuity as the limit of a thin boundary 
layer over which the properties of the tube vary. Specifically 
we assume 

{1 if r < To, 

(ri -r)/(ri -ro) ifro < r < ri, (21) 
if r-i < r , 

where r is the distance from the tube axis, and consider the 
limit ri vq. In practice we resolve ri — tq using 5 grid 
points and increase the resolution until the behaviour con- 
verges. In any case we choose tq and ri such that the total 
flux corresponds to a 2 Mm tube with field strength -Btubc- 
This completes the specification of the background. On 
this background we consider a plane parallel wave. The test 
cases considered here assume that the direction of propa- 
gation is perpendicular to the flux tube so that the prob- 
lem is essentially two-dimensional. Since we use different 
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Fig. 3 Scattered pressure perturbations (snapshot) from a 
horizontal 1 kG flux tube with axis x = z = 0. The radius 
of the tube is again 2 Mm. The bottom panel shows the exact 
result and the top panel the numerical result. 

schemes in the horizontal and vertical dkections we have 
checked the code using vertical and horizontal tubes. Fig- 
ures [T] and |2] show a comparison of the scattered pressure 
perturbations computed analytically and numerically for a 
vertical 1 kG flux tube. The numerical calculation used 200 
Fourier modes in both the x and y directions. Similar results 
were found for the horizontal tube (see Fig [3]). 

It is important to comment that the resolution (24 Mm 
over 200 pixels in each direction) is required to accurately 
capture the wave interaction with the discontinuities at the 
edge of the tube rather than the tube itself. For problems 
without discontinuities we expect convergence will occur at 
even lower resolutions. 

It is also worthwhile to show the results at a lower res- 
olution (here 100 pixels in each direction) with only three 
points resolving the boundary layer (ri — ro). In this case we 
have chosen a strong 3 kG magnetic field strength for which 
the density in the flux-tube has reduced by between 98% 
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Fig. 4 Scattered pressure perturbations (snapshot) from a 
strong, 3 kG, vertical flux tube with axis x ~ y — 0. The 
bottom panel shows the exact result, the top panel the nu- 
merical. The calculation was performed at a lower resolu- 
tion of 100 by 100 grid points with only 3 points across the 
tube boundary. This shows the type of artifacts which can 
be expected from using a too low resolution. 

and 99%. The results are shown in Figures |4] and |5] Mag- 
netic fields of this strength, and this level of evacuation, will 
presumably be important in calculations of wave propaga- 
tion through sunspots. Given the resolution, the agreement 
between the numerical and analytic results are reasonable, 
however artifacts near the now poorly resolved boundary 
layer are clearly visible. 

4.2 Background velocity field 

The above test involves several terms related to background 
inhomogeneities in the magnetic field, the density and the 
pressure fields. An obvious remaining class of perturbations 
involve the velocity field, and we have tested the code here 
by comparing against an analytic solution for uniform flow 
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Fig. 5 A cut at y = through the scattered pressure per- 
turbations shown in Figure |4] (vertical 3 kG flux tube, low 
resolution). The dashed line shows the exact solution and 
the solid line the numerical. The units are such that the in- 
coming pressure perturbations have an amplitude of 1 . 

in a cylinder embedded in a stationary uniform media. We 
assume that the density and temperature are the same inside 
and outside the cylinder, so that the only defining property 
is the velocity which we assume is directed along the axis of 
the cylinder The analytic solution was obtained in a similar 
manner to that in Gizon et al. (2006). 

An important property of the solution is that the scat- 
tered wave field vanishes for a plane parallel wave perpen- 
dicular to the cylinder so that all non-trivial solutions are 
three-dimensional. We have therefore chosen to compare 
the analytic and numerical solutions for the case where the 
cylinder and direction of wave propagation make a 45° an- 
gle. An additional complication is that the analytic solution 
no longer corresponds with an undisturbed plane wave en- 
countering an obstacle: rather part of the wave-front has pre- 
viously passed through the cylinder 

For the calculation we have taken the wave propagation 
is in the x direction with the cylinder running from the bot- 
tom left to the top right in the x-z plane. The radius of the 
cylinder is 2 Mm, the background sound speed is 10 km/s, 
and the flow in the cylinder is 1 km/s. Again there is a sharp 
jump between the properties inside the cylinder and those 
outside, and in principle we should use different resolu- 
tions to check the convergence of the solution. In practice 
we have used a grid with 200 grid points or modes in each 
direction, a similar resolution to that found suitable for the 
magnetic/pressure perturbations. 

The analytic and numerical solutions were found to be 
in reasonable agreement, as is shown in Figures. |6] and [T] 
Unlike the results shown in section 4.1, an exact agreement 
in this case cannot be expected since the numerical and ana- 
lytic problems are significantly different. The analytic solu- 
tion is a snapshot of a wave which has always been affected 
by the tube. In contrast the numerical solution evolves, with 
the scattering beginning at t = t^. What makes these two 
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solutions a useful test is that the two problems become iden- 
tical in the limit of a large box. If we kept increasing the 
size of our box the analytic and numerical solutions would 
presumably converge. Tests with smaller boxes (cubes of 
12 Mm) indeed showed a much worse match between the 
analytic and numerical results, so the agreement we are 
achieving with even 24 Mm boxes is reasonable. 

5 Discussion 

We have written and tested a code for studying wave prop- 
agation in an inhomogeneous magnetised solar atmosphere. 
There is a wide range of possible appUcations of the code, 
including understanding the observational signature of wave 
propagation through small-scale magnetic features, through 
granulation and supergranulation, and through the complex 
magnetic and velocity fields associated with sunspots and 
their penumbrae. Tests of the SLiM code so far have been 
successful. We have confidence that the code will be a use- 
ful tool to understand the signatures of wave propagation 
through various types of solar features. 
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Fig. 6 Scattered pressure perturbations (snapshot) from a 
1 km/s flow within a cylinder inclined at 45° to the the di- 
rection of propagation of the incoming wave. The flow is 
aligned with the axis of the cylinder, directed upwards and 
to the right. Shown is the x-z plane. The black lines indi- 
cate the boundaries of the cylinder. The bottom panel shows 
the analytic result, the top corresponds to the numerical so- 
lution. The initial condition for the numerical and analytic 
treatments are different, with the analytic treatment consid- 
ering a wave which has been in permanent contact with the 
cylinder. In the top panel, the lower left hand corner is least 
contaminated by the resulting artefacts. Note that because 
of the periodic nature of the box scattered waves appear in 
the top half of the box. The white line shows the cut used 
for Figure |7] 
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Fig. 7 A cut across the white Hnes shown in Figure|6] The 
solid curve shows the numerical solution, the dashed curve 
the analytic. 
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